Pfaffian pairing wave functions in electronic structure quantum Monte Carlo 
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We investigate the accuracy of trial wave function for quantum Monte Carlo based on pfafhan 
functional form with singlet and triplet pairing. Using a set of first row atoms and molecules we 
find that these wave functions provide very consistent and systematic behavior in recovering the 
correlation energies on the level of 95%. In order to get beyond this limit we explore the possibil- 
ities of multi-pfaffian pairing wave functions. We show that a small number of pfaffians recovers 
another large fraction of the missing correlation energy comparable to the larger-scale configuration 
interaction wave functions. We also find that pfaffians lead to substantial improvements in fermion 
nodes when compared to Hartree-Fock wave functions. 

PACS numbers: 02.70.Ss, 71.15.Nc, 71.15.-m 



Wave functions of interacting quantum systems such 
as electrons in matter are notoriously difficult to calcu- 
late despite decades of effort. One of the most productive 
many-body methods for electronic structure problems is 
quantum Monte Carlo (QMC) which employs stochastic 
approaches both for solving the stationary Schrodinger 
equation and for evaluation of expectation values Q,Q>l3- 
The key advantage of QMC is its capability to use ex- 
plicitly correlated wave functions which allow to study 
many-body effects beyond the reach of mean-field ap- 
proaches. The most important limit on QMC accuracy 
is the fixed-node approximation which is used to circum- 
vent the fermion sign problem 0, 0. Even with this 
approximation, the fixed-node QMC has been very suc- 
cessful for a host of real systems such as molecules, clus- 
ters and solids with up to hundreds of valence electrons 
and has provided an agreement with experiments (1-3%) 
for cohesive energies, band gaps, and other energy dif- 
ferences However, reaching beyond the fixed-node 
limit has proven to be challenging since fermion nodes 
(subset of position space where the wave function van- 
ishes) are complicated high-dimensional manifolds which 
are affected by correlation as well. Accuracy of QMC 
trial wave functions is therefore crucial for both funda- 
mental and also computational reasons. 

The state-of-the-art QMC calculations employ accu- 
rate Slater-Jastrow wave functions that can be written as 
vJr T = \g A exp[[/ corr ] where is the antisymmetric part 
while U corr describes the electron-electron and higher- 
order correlations. The antisymmetric part is either a 
single Hartree-Fock (HF) determinant of one-particle or- 
bitals or a multi-reference sum of excited determinants 
such as a limited Configuration Interaction (CI) expan- 
sion [lfjj . A natural generalization of one-particle orbital 
is a two-particle or pair orbital, sometimes called a gem- 
inal. In particular, the Bardeen-Cooper-Schrieffer (BCS) 
wave function, which is an antisymmetrized product of 



singlet pairs, has been recently used to calculate several 
atoms and molecules as well as superfluid Fermi gases 
[hi Hrl Hfl | . The results show significant improvements 
over the single-determinant HF wave functions. Never- 
theless, for partially spin-polarized systems the energy 
gains are less pronounced due to the lack of pairing cor- 
relations in the spin-polarized subspace [l^. The spin- 
polarized triplet pairing wave functions based on p faffians 
have been tried a few times on model systems |l ll Il2l Il8| . 

In this letter, we propose to describe systems of elec- 
trons by pfaffian wave functions with variational free- 
dom beyond HF and also BCS wavefunctions. The pfaf- 
fian allows us to incorporate pair orbitals for both singlet 
and triplet pairing channels together with unpaired one- 
particle orbitals into a single, compact wave function. 
These pfaffian wave functions are tested on atomic and 
molecular systems in variational and fixed-node diffusion 
Monte Carlo methods. The results show significant gains 
in correlation energy both for spin-polarized and unpo- 
larized cases. Furthermore, we explore the multi-pfaffian 
wave functions and we find that they recover a large frac- 
tion of the missing correlation energy while being much 
more compact than expansions in determinants. 

Let us consider 1,2, ...,N spin-up and N + 1,...,2N 
spin-down electrons in a singlet state with electron spa- 
tial coordinates denoted simply as R = (1, 2, 2iV). An 
antisymmetrized product of singlet pair orbitals <fr(i,j) — 
4>(j, i) is the BCS wave function 



BCS 



= A[<l>(i,j)] =dd#(i,j)] = det[*] 



(1) 



which is simply a determinant of N x N matrix. The 
BCS wave function is efficient for describing systems with 
single-band correlations such as Cooper pairs in con- 
ventional BCS superconductors where pairs form from 
one-particle states close to the Fermi level. For par- 
tially spin-polarized states one can augment the matrix 
by columns/rows of one-particle orbitals. However, the 
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spin-polarized subspace is then uncorrelated and for a 
fully spin-polarized system one ends up with the usual 
Hartree-Fock wave function. In order to correlate spin- 
polarized electrons it is necessary to generalize the wave 
function form and introduce effects of triplet pairing. 

For a system of 2N fully spin-polarized electrons the 
pairing wave function is formed as an antisymmetrizcd 
product of triplet pair orbitals x(*>j) — ~x{hi) and is 
given by [II El El 



A[ X {1, 2) X (3, 4)...] - pf [ X (i, j)] = pf M 
which defines a pfaffian of degree 27V, eg, for N = 2 



(2) 



pf Mhj)} = X(l, 2)x(3, 4)-%(l, 3) X (2, 4)+ X (l, 4) X (2, 3). 

(3) 

Any pfaffian of an odd degree vanishes, however, the 
pfaffian wave function can be easily generalized to an 
odd number of electrons by extending the pfaffian by 
a row/column of one-particle (unpaired) orbital. For 
example, for three spin-up electrons, replace the last 
row/column in the equations above as x(*?4) — > ip(i) 
and x(4,«) — > —<p(i). We note that the Hartree-Fock 
wave function is a special case of both BCS singlet and 
pfaffian triplet pairing wave functions as discussed |2(1 |. 

The square of the pfaffian is related to the determinant 
of a skew-symmetric matrix as 



{pf[x(*,j)]} 2 =dct[ X (z,j)]- 



(4) 



However, the QMC applications require also the knowl- 
edge of the wave function sign, eg, for enforcing the 
fixed-node restriction. Therefore, we have implemented 
a direct evaluation of pfaffian based on an 0(N 3 ) algo- 
rithm which is analogous to Gauss elimination for deter- 
minants. Note that pfaffians can be expanded in pfaffian 
minors and by exploring Cayley's results [9j one can cal- 
culate the pfaffian and its updates for electron moves in 
computer time similar to calculation of determinants |20| . 

Let us now consider a partially spin-polarized system 
with unpaired electrons. Remarkably, the pfaffian form 
can accommodate both singlet and triplet pairs as well 
as one-particle unpaired orbitals into a single, compact 
wave function. The singlct/triplct/unpaired (STU) or- 
bital pfaffian wave function is given by 







<j>Tl 






*STC7 = Pf 








(5) 
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where the bold symbols are block matrices/vectors of cor- 
responding orbitals and T denotes transposition. For 
a spin-restricted STU wave function the pair and one- 
particle orbitals of spin-up and -down channels would be 
identical. 

The pfaffian wave functions were used in QMC cal- 
culations by variational and fixed-node diffusion Monte 



Carlo (VMC and DMC) methods [a |g. The VMC 
trial/ variational wave function is a product of an anti- 
symmetric part ^ a times a Jastrow correlation factor 

^VMC(R) = ^A(R)e-Kp[ U corr({rij},{ r il},{ r ij})} (6) 

where U corr depends on electron-electron, electron-ion 
and electron-electron-ion combinations of distances 0, 
llij with maximum of 22 variational parameters. For 
the antisymmetric part we have used ^a — ^hf and 
^a = ^stu as well as some tests with ^a = ^bcs to 
compare with recent results 0, 0] . The pair orbitals 
were expanded in products of one-particle orbital basis 

EF 



as 



(i,j),x(i,j) 



k,l 



Ckitpk(i)<Pi(j)- 



(7) 



The coefficients are symmetric (cki = cu~) for the singlet 
4>{i,j), and antisymmetric {cm — — c;fc) for the triplet 
x{hj) functions. The expansions include both occupied 
and unoccupied (virtual) one-particle orbitals. The one- 
particle atomic and molecular orbitals used in expan- 
sions, which we tested, were either Hartree-Fock orbitals 
or natural orbitals |l(j from CI correlated calculations. 
Typically, we used about 10 virtual orbitals and the nat- 
ural orbitals produced better and more systematic results 
than the HF ones. The pair orbital expansion coefficients 
were then optimized in VMC by minimizations of energy 
using recently published method |l8| . 

We have applied these developments to several first row 
atoms and dimers (Fig. Except for the Be atom, we 
used pseudopotentials to eliminate the atomic cores E3 
while the previous calculations with BCS wave functions 
were done with all electrons |l5j. Nevertheless, our BCS 
wave functions produced percentages of correlation ener- 
gies rather close to the ones obtained with all electrons 
Perhaps the most striking result is a systematic 
percentage of recovered correlation energy (94-97%) for 
systems heavier than Be (see Fig. ^| . The triplet contri- 
bution for these single pfaffian STU wave functions are 
small, with the only exception being the nitrogen atom, 
where we see a gain of additional I % in correlation en- 
ergy (Tabic P) when compared to the trial wave function 
without triplet pairs. We believe that this is due to the 
fact that the nitrogen atom is quartet and therefore has 
the highest spin polarization from all studied cases. 

Our results show that the single pfaffian form is capa- 
ble of capturing single-band near-degeneracies and mix- 
ing of excited states for both spin-polarized and unpolar- 
ized systems. Considering multi-dctcrminantal expan- 
sions, such as the CI method, the overall trade-off be- 
tween accuracy and computational cost seems to be in 
favor of more compact and physically-based pfaffian wave 
functions. A similar opinion was also expressed by Sorella 
and coworkers 151 . 
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FIG. 1: Correlation energies obtained by QMC methods with 
the different trial wave functions: VMC and fixed-node DMC 
with HF nodes (HF) and STU pfaffian nodes (PF). The lower 
plot shows the fixed-node DMC correlation energy gains over 
HF nodes for BCS and STU pfaffian wave functions. The sta- 
tistical error bars are of the symbol sizes or smaller. Except 
for the Be atom all the calculations used the same pseudopo- 
tentials p^|. 



TABLE I: Total energies (a.u.) for N atom and dimer with 
amounts of correlation energy recovered, in VMC and DMC 
methods with wave functions as discussed in the text. 



WF 



N 



N 2 



Ecor?' 



HF 

VMC/HF 

VMC/BCS 

VMC/STU 

DMC/HF 

DMC/BCS 

DMC/STU 

Exact /est. 



-9.628915 
-9.7375(1) 
-9.7427(3) 
-9.7433(1) 
-9.7496(2) 
-9.7536(2) 
-9.7551(2) 
-9.759215 





83.3(1) 
87.3(2) 
87.8(1) 
92.6(2) 
95.7(2) 
96.8(1) 
100 



-19.44946 
-19.7958(5) 
-19.8179(6) 

-19.821(1) 
-19.8521(3) 
-19.8605(6) 
-19.8607(4) 

-19.88196 





80.1(1) 
85.2(1) 
86.0(2) 
93.1(1) 
95.1(1) 
95.2(1) 
100 



X3[2s,2p x ]. We can therefore include all three pfaffians 
into our ^mpf and further optimize independently all 
the pairing functions in VMC on the space of occupied 
and virtual orbitals. This construction allows us to in- 
corporate the excitations which are not present in the 
single pfaffian STU wavefunction based only on a single 
$ and a single x pairing functions. In the leading order 
the resulting MPF wave function then corresponds to the 
CI wave function with singles and doubles with the same 
active orbital space. The disadvantage of this approach 
is that we perform the VMC optimizations of M 2 pairing 
coefficients for each pfaffian given that M is the total size 
of our orbital basis. However, this can be improved by a 
factor of M if we rotate the one-particle orbitals to make 
the pairing functions diagonal. The total number of pfaf- 
fians in the expansion is then subject to required sym- 
metry of state and desired accuracy. In the most general 
case it would be proportional to the number of distinct 
pairs. For large systems one can restrict the number of 
pairs by including only the ones which have significant 
contributions, eg, pairs formed from states close to the 
Fermi level, or considering pairs only within a given band 
or a sub-band. 



TABLE II: Percentages of correlation energies recovered for 
C, N and O atoms by VMC and DMC methods with wave 
functions as discussed in the text. Corresponding number of 
pfaffians/determinants n for each wave function is also shown. 
The estimated exact correlation energies for C,N,0 are 0.1031, 
0.1303, 0.1937 a.u. 



WF 


n C 


n 


N 


n 


O 


VMC(MPF) 


3 92.3(1) 


5 


90.6(1) 


11 


93.6(2) 


VMC(CI) 


98 89.7(4) 


85 


91.9(2) 


136 


89.7(4) 


DMC(MPF) 


3 98.9(2) 


5 


98.4(1) 


11 


97.5(1) 


DMC(CI) 


98 99.3(3) 


85 


98.9(2) 


136 


98.4(2) 



In order to test the limits of pfaffian functional form, 
we propose a simple extension: a multiple pfaffian (MPF) 
wave function having a form 

^mpf = pfbd, <t>\, ft] + pf[X2, <f>2, <P2] + ■ ■ ■ (8) 

so that in Eq. HJwe have ^>a = ^mpf- In actual calcu- 
lations we start with all pairing functions such that each 
pfaffian is equal to the HF wave function, pf <fii, <pi] = 
^hf- The pairing orbitals [see Eq. Q] are expanded in 
the basis of HF/natural occupied orbitals, eg, for the car- 
bon atom we have 2s, 2p x and 2p y . The choice of singlet 
0i(l, 2) = 2s(l)2s(2) = <M2s,2s] and triplet xi(l»2) = 
2p x (l)2p y (2) - 2p y (l)2p x (2) = xi[2px,2p v ] pairing or- 
bitals then gives pf[xi,</>i] = ^hf[2s^ 7 2p x , 2p y \. How- 
ever, one can construct the equivalent combinations 
of pairs as: 4> 2 [2s,2p x ), X2[2s,2p y ] and (f> 3 [2s,2p y ], 



The results (Table [Fiji show that our MPF wave func- 
tions are able to recover close to 99% of correlation en- 
ergy. Furthermore, comparison with the CI results shows 
that it is possible to obtain similar quality of wave func- 
tions with corresponding improvements of the fermion 
nodes at much smaller calculational cost. This is another 
indication that the inclusion of singlet and triplet pairing 
enables us to treat the correlation in both spin-polarized 
and -unpolarized channels in a consistent manner. 

The quality of fermion nodes is crucial for accurate 
energies in the fixed- node DMC. The fermion node man- 
ifold is defined by an implicit equation ^(R) = 0. For 
N electrons the manifold has (3N — 1) dimensions and 
divides the configuration space into compact nodal cells. 
The HF spin-polarized ground states typically show two 
nodal cells, while HF unpolarized or partially polarized 
states have 2x2 = 4 nodal cells since the wave function 
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only the lower bound on the recovered correlation en- 
ergies since we were focused on systematic trends using 
new functional forms. For example, a simultaneous reop- 
timization of one-particle orbitals used in representation 
of pairs could improve the results further. Finally, the 
pfaffian wave functions exhibit qualitative improvements 
in fermion nodes and eliminate a significant portion of 
the Hartree-Fock node errors. 

We gratefully acknowledge the support by ONR- 
N00014-01-1-0408, and NSF DMR-0121361, DMR- 
0121361, and EAR-0530110 grants and the computer 
time allocations at NCSA, PSC and SDSC facilities. 



FIG. 2: A 3D cut through the fermion node hypersurface 
of oxygen atom obtained by scanning the wave function with 
a pair of spin-up and -down of electrons, both sitting at the 
scanning point. The remaining electrons are fixed at a given 
VMC snaphsot positions (small green spheres). Nucleus is 
depicted in the center of the cube by the black sphere. The 
three colors show nodes of: Hartree-Fock (red/dark gray); 
STU pfaffian nodes (orange/medium gray); and the nodes of 
the CI wave function (yellow/light gray). The CI node is 
very close to the exact one (see text). The HF node clearly 
divides the space into four nodal cells while pfaffian and CI 
wave functions partitioning leads to the minimal number of 
two nodal cells. 

is a product of spin- up and -down determinants 0, 0, @ ■ 
Changes in the HF nodal structures from inclusion of cor- 
relation were recently investigated by D. Bressanini et al. 
0, H| • Here we observe that the pairing correlations have 
an important effect on the nodal structure as well. A 
changes in the nodal manifold topology is illustrated in 
Fig. [21 on the example of oxygen atom. As expected, the 
HF nodes show four nodal cells while the pfaffian " opens" 
these artificial compartments and changes the topology 
to the minimal number of two nodal cells, similar to the 
effect of correlation observed for Be atom Q ■ A compar- 
ison of pfaffian nodes with very accurate CI nodes shows 
that both are qualitatively similar. Further results on 
the properties of nodes and their improvement are given 
elsewhere 

In conclusion, we have proposed pfaffians with sin- 
glet pair, triplet pair and unpaired orbitals as varia- 
tionally rich and compact wave functions which offer 
significant and systematic improvements over commonly 
used Slater determinant-based wave functions. We have 
demonstrated that these pfaffian pairing wave functions 
are able to capture a large fraction of missing correlation 
energy with consistent treatment of correlation for both 
spin-polarized and unpolarized pairs. We have explored 
also multi-pfaffian wave functions which enable us to ob- 
tain more correlation while keeping the wave functions 
compact. Our pfaffian wave function results represent 
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